#!/usr/bin/env python
# -*- coding: utf-8 -*-

# Copyright (c) 2011 - 2013 Stefano Mazzucco <stefano -at- curso.re>
# All rights reserved.
#
# This file is part of Crystal Ball Plus.
#
# Crystal Ball Plus is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# Crystal Ball Plus is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with Crystal Ball Plus.  If not, see <http://www.gnu.org/licenses/>.

"""Cristallographic classes.

"""

import numpy as np
import os

import spacegroups as sgs
from fileio import nl, newline, comment, special, overwrite, newname
from operations import *


def _is_deg(unit):
    """Return True if unit is 'deg', 'degree', 'degrees',
    return False if unit is 'rad', 'radian', 'radians',
    raise ValueError otherwise.

    Case insensitive.

    """
    deg = None
    if unit.lower() in ('deg', 'degree', 'degrees',):
        deg = True
    elif unit.lower() in ('rad', 'radian', 'radians',):
        deg = False
    else:
        msg = "'%s' is neither 'degrees' nor 'radians'" % unit
        raise ValueError(msg)
    return deg

class MillerPlane(object):
    """Class that defines the Miller indices of a crystallographic plane.

    A MillerPlane instance must be declared as e.g.

       >>> m = MillerPlane(plane) # 'plane' is a 3-iterable

    """
    def __init__(self, hkl):
        err = ''
        hkl_dict = {0: 'h', 1: 'k', 2: 'L'}
        for i in xrange(len(hkl)):
            if int(hkl[i]) != hkl[i]:
                err += '\nIndex %s = %s is not integer' % (hkl_dict[i], hkl[i])
        if err:
            raise TypeError(err)
        self.hkl = np.asarray(hkl, dtype='int64')
        if len(self.hkl.shape) != 1 or (self.hkl.shape[0] != 3):
            err = '%s indices instead of 3' % len(self.hkl.shape)
            raise AttributeError(err)

    def __iter__(self):
        return self.hkl.__iter__()

    def h(self):
        """Return the value of the first Miller index 'h'

        """
        return self.hkl[0]
    def k(self):
        """Return the value of the second Miller index 'k'

        """
        return self.hkl[1]
    def L(self):
        """Return the value of the first Miller index 'L'

        .. note::
           In compliance with PEP8, an uppercase letter
           is used instead of a lowercase to avoid confusion between
           the letter 'l' and the number '1'

        """
        return self.hkl[2]

    def angle(self, hkl, uc, deg=True):
        """Return the angle with the Miller indices defined by hkl
        in the unit cell definde by uc.

        *Parameters*

         hkl : instance of MillerPlane or array of 3 elements

        uc : array of 6 elements
            Unit cell: (a, b, c, alpha, beta, gamma)

        deg : boolean
              whether the angles in the unit cell are degrees (True, default)
              or radians (False)

        *Returns*

        theta : float
                angle with plane (hkl)
                if deg=True (default) the units will be degrees,
                if deg=False the units will be radians

        """
        if hkl.__class__.__name__ is 'MillerPlane':
            plane = hkl.hkl
        else:
            plane = hkl
        mt = reciprocal_metric_tensor(uc, deg)
        return theta(self.hkl, plane, mt, deg)

    def spacing(self, uc, deg=True):
        """Return the value of the d-spacing corresponding in the unit cell
        defined by uc.

        *Parameters*

        uc : array of 6 elements
            Unit cell: (a, b, c, alpha, beta, gamma)

        deg : boolean
              whether the angles in the unit cell are degrees (True, default)
              or radians (False)

        *Returns*

        d : float
            d-spacing of the current MillerPlane instance

        """
        return d_spacing(self.hkl, uc, deg)

    def zone_axis(self, hkl, deg=True):
        """Return the zone axis [u, v, w] with the crystallographic planes
        determined by the Miller indices hkl.

        *Parameters*


        hkl1 :  instance of MillerPlane or array of 3 elements

        deg : boolean (optional)
              Whether the angles in the unit cell are degrees (True, default)
              or radians (False)

        *Returns*


        uvw : array of 3 elements
              zone axis (u, v, w)

        """
        if hkl.__class__.__name__ is 'MillerPlane':
            plane = hkl.hkl
        else:
            plane = hkl
        return zone_axis(self.hkl, plane, deg)

    def family(self, sg,  uc, deg=True):
        """Return the set of equivalent directions (family) in the crystal
        system defined by the space group 'sg' and unit cell 'uc'.

        *Parameters*

        sg : instance of SpaceGroup

        uc : array of 6 elements
             Unit cell: (a, b, c, alpha, beta, gamma)

        deg : boolean (optional)
              Whether the angles in the unit cell are degrees (True, default)
              or radians (False)

        *Returns*

        family : list of MillerPlane

        *See also*

        family, Diffractogram.family

        """
        mlist = family(self.hkl, sg, space='r', uc=uc, deg=deg)
        for i, m in enumerate(mlist):
            mlist[i] = MillerPlane(m)
        return mlist

    def __eq__(self, obj):
        if obj.__class__.__name__ is 'MillerPlane':
            return self.hkl.__eq__(obj.hkl)
        else:
            return self.hkl.__eq__(obj)

    def __repr__(self):
        out = self.__class__.__name__ + '('
        for i in self.hkl:
            if i < 0:
                out += repr(i) + ', '
            else:
                out += ' ' + repr(i) + ', '
        out = out.rstrip(' ,')
        out += ')'
        return out

    def __str__(self):
        out = self.__class__.__name__ + '('
        for i in self.hkl:
            if i < 0:
                out += repr(i) + ' '
            else:
                out += ' ' + repr(i) + ' '
        out = out.rstrip()
        out += ')'
        return out

class Diffractogram(object):
    """A generic diffractogram.

    Be particularly careful about the units of the angles and spacings!
    Usually, they are given in degrees and Angstroms respectively.

    The keyword 'fp' can be either a file or a dictionary with the
    following keywords::

       {'chemname' :  string
       'structure' : string
       'sg_name' : string
       'sg_number' : integer
       'crystsys' : string
       'unitcell' : array of 6 elements
       'uc_len_units' : string
       'uc_ang_units' : string}

    Example

       >>> d = Diffractogram('/path/to/some/file.dfg')

    """
    ext = '.dfg'                # diffractogram file extension
    def __init__(self, fp='', sep=None):
        self.filename = ''

        self.sep = sep          # separator

        self.structure = ''         # e.g. Fe3C
        self.chemname = ''          # e.g. Iron Carbide
        self.sg_name = ''           # e.g. PNMA
        self.sg_number = None       # e.g. 225
        self.sg = None              # instance of class SpaceGroup
        self.crystsys = ''          # e.g. tetragonal
        self.unitcell = np.zeros((6,)) # 6-tuple: (a, b, c, alpha, beta, gamma)
        self.uc_len_units = ''
        self.uc_ang_units = ''

        self.dmt = None         # direct metric tensor
        self.rmt = None         # reciprocal metric tensor

        self.hkl = None         # Miller indices

        self.d_spacings = None
        self.d_spacings_units = ''       # usually, it's Angstroms

        self.angles_to_x = None      # angles to horizontal axis
        self.angles_to_x_units = ''   # usually, it's degrees

        self.intensities = None
        self.intensities_units = ''

        self.image = None

        self.info = {
            'file name' : self.filename,
            'chemical name' : self.chemname,
            'structure' : self.structure,
            'space group' : self.sg_name,
            'crystal system' : self.crystsys,
            'unit cell' : self.unitcell,
            }
        if fp:
            if hasattr(fp, 'has_key'):
                for key, value in fp.iteritems():
                    # setattr could be used as well.
                    if key == 'chemname' :  self.chemname = value
                    if key == 'structure' : self.structure = value
                    if key == 'sg_name' : self.sg_name = value
                    if key == 'sg_number' : self.sg_number = value
                    if key == 'crystsys' : self.crystsys = value
                    if key == 'unitcell' : self.unitcell = value
                    if key == 'uc_len_units' : self.uc_len_units = value
                    if key == 'uc_ang_units' : self.uc_ang_units = value
            elif os.path.exists(fp):
                self.filename = os.path.basename(fp)
                with open(fp, 'r') as f: # must be a text file
                    if os.path.splitext(fp)[1].lower() != Diffractogram.ext:
                        msg = '\nIt looks like your file has'
                        msg += (' a wrong extension ("%s"'
                                % os.path.splitext(fp)[1])
                        msg += ' insetad of "%s").\n' % Diffractogram.ext
                        msg += 'The data might not be read properly.\n'
                        print(msg)
                    print("Reading diffractogram from file '%s'" % self.filename)
                    self.read_file(f)
            else:
                err = 'File %s does not exist.' % fp
                raise IOError(err)
                return
        if not np.any(self.d_spacings):
            # if there are no d-spacing, generate some
            self.gen_d_spacings(start=0, stop=2)
        if np.any(self.intensities):
            self.sort_by_int()
        msg = ''
        if self.sg_name:
            try:
                # trigonal hexagonal/rombohedral structure ?
                # http://cci.lbl.gov/~rwgk/my_papers/CCN_2011_01_H3_H32.pdf
                if (self.sg_name[-1].lower() == 'r' or self.crystsys.lower().find('rombohedral') != -1) and not ( np.all(self.unitcell[3:] == (90., 90., 120.)) or np.all(self.unitcell[3:] == (np.pi/2., np.pi/2, np.pi*2/3))):
                    if self.sg_number:
                        self.sg_number = 1000 + self.sg_number
                    self.sg = sgs.get_spacegroup(self.sg_name[:-1])
                else:
                    self.sg = sgs.get_spacegroup(self.sg_name)
            except ValueError, msg:
                pass
        if msg and self.sg_number:
            msg = ''
            try:
                self.sg = sgs.get_spacegroup(self.sg_number)
            except ValueError, msg:
                pass
        elif self.sg_number:
            try:
                self.sg = sgs.get_spacegroup(self.sg_number)
            except ValueError, msg:
                pass
        if msg:
            print('SpaceGroup instance not created (%s)' % msg)
        if self.sg is not None:
            self.sg_number = self.sg.number
            self.sg_name = self.sg.short_name
            self.crystsys = self.sg.crystal_system.lower()
        self.update_info()
        self.update_dmt()
        self.update_rmt()

    def update_dmt(self):
        """Update the direct metric tensor attribute dmt.

        """
        if np.any(self.unitcell):
            try:
                deg = _is_deg(self.uc_ang_units)
            except ValueError, msg:
                print(msg)
                print('Assuming degrees')
                deg = True
            self.dmt = direct_metric_tensor(self.unitcell, deg)

    def update_rmt(self):
        """Update the reciprocal metric tensor attribute rmt.

        """
        if np.any(self.unitcell):
            try:
                deg = _is_deg(self.uc_ang_units)
            except ValueError, msg:
                print(msg)
                print('Assuming degrees')
                deg = True
            self.rmt = reciprocal_metric_tensor(self.unitcell, deg)

    def update_info(self):
        """Update the info attribute.

        """
        self.info = {
            'file name' : self.filename,
            'chemical name' : self.chemname,
            'structure' : self.structure,
            'space group' : self.sg_name,
            'crystal system' : self.crystsys,
            'unit cell' : self.unitcell,
            }

    def family(self, v, space='d', rnd=False):
        """Return the family of directions/planes equivalent to 'v'.

        *Parameters*

        v : array of 3 elements or instance of MillerPlane

        space : string (optional)
                space can be 'd' (direct, default) or 'r' ('reciprocal')
                if v is an instance of MillerPlane, space is automatically
                set to 'r'

        rnd : boolean (optional)
              whether to round the result using numpy.around (rnd=False,
              default) or not (rnd=True)
              NOTE: if space='r', rnd is set to 'True'

        *Returns*

        family : list of arrays

        *See also*

        family, MillerPlane.family

        """
        if self.sg is None:
            raise ValueError("space group 'sg' is %s" % self.sg)
        if not np.all(self.unitcell):
            raise ValueError("unit cell 'unitcell' is %s" % self.unitcell)
        try:
            deg = _is_deg(self.uc_ang_units)
        except ValueError, msg:
            print(msg)
            print('Assuming degrees')
            deg = True
        if v.__class__.__name__ is 'MillerPlane':
            return v.family(self.sg, self.unitcell, deg=deg)

        return family(v, self.sg, space=space, uc=self.unitcell,
                      deg=deg, rnd=rnd)

    def stored_spacing(self, *hkl):
        """Return the value of the stored d-spacing (if any) corresponding
        to the given Miller indices hkl.

        *Parameters*

        *hkl : 3 integers or array of 3 elements or instance of MillerPlane

        Examples:
           >>> d = Diffractogram('some_file')
           >>> d.stored_spacing(0, 0, 2)
           >>> a = (1, -1, 0)
           >>> d.stored_spacing(a)
           >>> b = MillerPlane([-1, 0, -1])
           >>> d.stored_spacing(b)

        """
        if len(hkl) == 1:
            if hkl[0].__class__.__name__ is 'MillerPlane':
                hkl = tuple(hkl[0].hkl)
            elif hasattr(hkl[0], '__iter__'):
                hkl = tuple(hkl[0])
        elif len(hkl) != 3:
            print('Invalid parameters.')
            return
        rows = []
        for i in xrange(self.hkl.shape[0]):
            rows.append(tuple(self.hkl[i].hkl))
        try:
            idx = rows.index(hkl)
            return self.d_spacings[idx]
        except:
            print('stored d-spacing corresponding to %s not found.' % (hkl,))
            return None

    def spacing(self, *hkl):
        """Return the value of the d-spacing corresponding to the given Miller
        indices hkl.

        *Parameters*

        *hkl : 3 integers or array of 3 elements or instance of MillerPlane

        Examples:

           >>> d = Diffractogram('some_file')
           >>> d.spacing(0, 0, 2)
           >>> a = (1, -1, 0)
           >>> d.spacing(a)
           >>> b = MillerPlane([-1, 0, -1])
           >>> d.spacing(b)

        """
        if not np.all(self.unitcell):
            print('Unit cell is not valid: %s' % self.unitcell)
            print('Nothing to do')
            return
        if len(hkl) == 1:
            if hkl[0].__class__.__name__ is 'MillerPlane':
                hkl = tuple(hkl[0].hkl)
            elif hasattr(hkl[0], '__iter__'):
                hkl = tuple(hkl[0])
        elif len(hkl) != 3:
            print('Invalid parameters.')
            return
        try:
            deg = _is_deg(self.uc_ang_units)
        except ValueError, msg:
            print(msg)
            print('Assuming degrees')
            deg = True
        return d_spacing((hkl), self.unitcell, deg)

    def angle(self, hkl_a, hkl_b):
        """Return the angle with the Miller indices defined by hkl_a
        and hkl_b in the unit cell definde by uc.

        *Parameters*

        hkl_a : instance of MillerPlane or array of 3 elements

        hkl_b : instance of MillerPlane or array of 3 elements

        *Returns*

        theta : float
                angle between plane hkl_a and hkl_b

        """
        if not np.all(self.unitcell):
            print('Unit cell is not valid: %s' % self.unitcell)
            print('Nothing to do')
            return
        if hkl_a.__class__.__name__ is 'MillerPlane':
            plane_a = hkl_a.hkl
        else:
            if len(hkl_a) != 3:
                print("Invalid parameters '%s'." % hkl_a)
                return
            plane_a = hkl_a

        if hkl_b.__class__.__name__ is 'MillerPlane':
            plane_b = hkl_b.hkl
        else:
            if len(hkl_b) != 3:
                print("Invalid parameters '%s'." % hkl_b)
                return
            plane_b = hkl_b
        try:
            deg = _is_deg(self.uc_ang_units)
        except ValueError, msg:
            print(msg)
            print('Assuming degrees')
            deg = True
        mt = reciprocal_metric_tensor(self.unitcell, deg)
        return theta(plane_a, plane_b, mt, deg)

    def gen_d_spacings(self, start=-5, stop=5, step=1):
        """Generate the d_spacings from the Miller planes and update the
        attributes 'd_spacings' and 'hkl'

        .. warning::
           Existing information about the intensity will not be valid
           anymore! Thus, the attribute `intensities` will be set to 'None'.

        By default, generate the d_spacings from all the valid planes from
        (-5, -5, -5) to (5, 5, 5).

        *Parameters*

        start : integer (optional)
                starting Miller index, default=0

        stop : integer (optional)
               ending Miller index, default=5

        step : integer (optional)
               step to advance to the next index

        """
        # TODO: exclude equivalent HKL depending on the crystal system
        # TODO: calculate intensities
        if not np.all(self.unitcell):
            # Nothing to do.
            return
        try:
            deg = _is_deg(self.uc_ang_units)
        except ValueError, msg:
            print(msg)
            print('Assuming degrees')
            deg = True
        newsp = []
        newhkl = []
        for h in xrange(start, stop + 1, step):
            for k in xrange(start, stop + 1, step):
                for L in xrange(start, stop, step):
                    if (h, k, L) != (0, 0, 0):
                        newsp.append(d_spacing((h, k, L), self.unitcell, deg))
                        newhkl.append(MillerPlane((h, k, L)))
        self.d_spacings = np.asarray(newsp, dtype='f4')
        self.hkl = np.asarray(newhkl, dtype='object')
        self.intensities = None

    def sort_by_int(self):
        """Sort the d-spacings and Miller indices by descending intensity.

        """
        if np.any(self.intensities):
            args = np.argsort(self.intensities)
            self.intensities = self.intensities[args]
            self.intensities = self.intensities[::-1]
            if np.any(self.d_spacings):
                self.d_spacings = self.d_spacings[args]
                self.d_spacings = self.d_spacings[::-1]
            if np.any(self.hkl):
                self.hkl = self.hkl[args]
                self.hkl = self.hkl[::-1]
            if np.any(self.angles_to_x):
                self.angles_to_x = self.angles_to_x[args]
                self.angles_to_x = self.angles_to_x[::-1]
        else:
            # Intensities are not present, cannot sort by intensity.
            return

    def read_file(self, fp):
        """Reads diffractogram data from file object 'fp' (must be text).
        The data are separated by separator 'sep' (defaults to whitespace).

        The file should have a (**case insensitive**) header of the form::

            ## d (d_unit) angle_to_x (angle_unit) hkl intensity

        At least one quantity among d, angle_to_x, hkl, and intensity
        is expected, the order does not matter. Units may be omitted.

        The valid keywords for the d-spacings column are:

        * d
        * d-spacing
        * d-spacings
        * d_spacing
        * d_spacings

        The valid keywords for the Miller planes column are:

        * hkl
        * miller

        .. note::

            In alternative to hkl, the header can have three separate columns
            'h', 'k', and 'l'.

        The valid keywords for the angles column are:

        * angle_to_x
        * angle-to-x
        * angles_to_x
        * angles-to-x
        * theta_to_x
        * theta-to-x

        The valid keywords for the intensities column are:

        * int
        * intensity
        * intensities

        Blank lines and comments (lines preceded by '#') are ignored.

        The following special comments (lines preceded by '##') are saved
        in memory::

            ## structure: string (e.g. Fe2 O3)

            ## chem name: string (e.g. Iron Oxide)

            ## space group: string (integer) (e.g. FM3-M (225))
               [values separated by 'sep']

            ## crystal system: string (e.g. orthorombic)

            ## unit cell: a b c alpha beta gamma (length_units) (angle_units)
               [values separated by 'sep']

        Class attributes that are created (if data is available):

        * d_spacings : array of float

        * d_spacings_units : string, usually 'A' or 'nm'

        * hkl : array of N instances of MillerPlane class

        * hkl_unit : string (useless, but you can specify it if you want)

        * angles_to_x : array of float (angles relative to horizontal axis)

        * angles_to_x_units : string, usually 'deg' or 'rad'

        * intensities : array of float

        * intensities_units : string, usually 'a.u.'

        * unitcell : 6-array of float (a, b, c, alpha, beta, gamma)

        * uc_len_units : string (e.g. 'A')

        * uc_ang_units : string (e.g. 'deg')

        * structure : string (e.g. 'H2O')

        * spacegroup : string (e.g. 'PNMA')

        """
        sep = self.sep
        if not fp.closed:
            # newline = ('\n', '\r\n')
            # comment = '#'
            # special = '##'

            struc = 'structure:'
            cn = 'chem name:'
            sg = 'space group:'
            cs = 'crystal system:'
            uc = 'unit cell:'

            d_s = ('d', 'd-spacing', 'd-spacings', 'd_spacing', 'd_spacings')
            hkl_s = ('hkl', 'miller')
            angle_s = ('angle_to_x', 'angle-to-x',
                       'angles_to_x', 'angles-to-x',
                       'theta_to_x', 'theta-to-x')
            int_s = ('int', 'intensity', 'intensities')
            h_s = ('h', )
            k_s = ('k', )
            L_s = ('l', )

            d_spacings = []
            hkl = []
            angles = []
            intensities = []

            d_col, hkl_col, angle_col, int_col = (None, None, None, None)
            h_col, k_col, L_col = (None, None, None)
            for line in fp.readlines():
                line = line.strip()
                if line[:2] == special:
                    line = line.lstrip(special).strip()
                    if struc in line:
                        self.structure = line[len(struc):].strip()
                    elif cn in line:
                        self.chemname = line[len(cn):].strip()
                    elif sg in line:
                        spacegroup = line[len(sg):].strip()
                        spacegroup = spacegroup.split(sep)
                        self.sg_name = spacegroup[0]
                        if len(spacegroup) == 2:
                            try:
                                self.sg_number = \
                                    int(spacegroup[1].strip('()'))
                            except ValueError, msg:
                                print('Space group number not set (%s)' %
                                      msg)
                    elif cs in line:
                        self.crystsys = line[len(cs):].strip()
                    elif uc in line:
                        unitcell = line.lstrip(uc)
                        unitcell = unitcell.split(sep)
                        if len(unitcell) == 6:
                            unitcell = np.asarray(unitcell, dtype='float')
                        elif len(unitcell) == 8:
                            self.uc_len_units = unitcell[6].strip('()')
                            self.uc_ang_units = unitcell[7].strip('()')
                            unitcell = np.asarray(unitcell[:6],
                                                  dtype='float')
                        else:
                            msg = 'Unit cell MUST have 6 elements'
                            msg += ' and, optionally, two unit strings'
                            raise IOError(msg)

                        self.unitcell = unitcell
                    else:
                        line = line.split(sep)
                        for word in line:
                            if word[0] == '(':
                                i = line.index(word)
                                line[i-1] = line[i-1] + line[i]
                                line.pop(i)
                        for word in line:
                            p = word.find('(')
                            if p == -1:
                                p = None
                            for s in d_s:
                                if s in word[:p].lower():
                                    d_col = line.index(word)
                                    if p:
                                        d_unit = word[p:]
                                        d_unit = d_unit.strip('()')
                                        self.d_spacings_units = d_unit
                                    break
                            for s in hkl_s:
                                if s in word[:p].lower():
                                    hkl_col = line.index(word)
                                    if p:
                                        hkl_unit = word[p:]
                                        hkl_unit = hkl_unit.strip('()')
                                        self.hkl_units = hkl_unit
                                    break
                            for s in h_s:
                                if s == word[:p].lower():
                                    h_col = line.index(word)
                                    break
                            for s in k_s:
                                if s == word[:p].lower():
                                    k_col = line.index(word)
                                    break
                            for s in L_s:
                                if s == word[:p].lower():
                                    L_col = line.index(word)
                                    break
                            for s in angle_s:
                                if s in word[:p].lower():
                                    angle_col = line.index(word)
                                    if p:
                                        angle_unit = word[p:]
                                        angle_unit = angle_unit.strip('()')
                                        self.angles_to_x_units = angle_unit
                                    break
                            for s in int_s:
                                if s in word[:p].lower():
                                    int_col = line.index(word)
                                    if p:
                                        int_unit = word[p:]
                                        int_unit = int_unit.strip('()')
                                        self.int_units = int_unit
                                    break
                else:
                    comment_idx =  line.find(comment)
                    if comment_idx >= 0 :
                        line = line[:comment_idx].strip()
                    if line:
                        line = line.split(sep)
                        if d_col is not None and line[d_col]:
                            d_spacings.append(line[d_col])
                        if hkl_col != None:
                            plane = line[hkl_col]
                            has_hkl = True
                        elif h_col != None and k_col != None \
                                and L_col != None:
                            plane = [line[h_col], line[k_col], line[L_col]]
                            has_hkl = True
                        else:
                            has_hkl = False
                        if has_hkl:
                            plane_3 = []
                            i = 0
                            while i < len(plane):
                                if plane[i] == '-':
                                    try:
                                        n = int(plane[i:i+2])
                                    except ValueError:
                                        n = None
                                    i = i + 2
                                else:
                                    try:
                                        n = int(plane[i])
                                    except ValueError:
                                        n = None
                                    i = i + 1
                                plane_3.append(n)
                            try:
                                plane_3 = MillerPlane(plane_3)
                            except TypeError:
                                msg = 'Skipping invalid hkl plane'
                                msg += ' from line %s' % line
                                print(msg)
                                plane_3 = '_invalid_'
                            hkl.append(plane_3)
                        if angle_col is not None and line[angle_col]:
                            angles.append(line[angle_col])
                        if int_col is not None and line[int_col]:
                            intensities.append(line[int_col])
            self.d_spacings = np.asarray(d_spacings, dtype='float')
            self.hkl = np.asarray(hkl)
            self.angles_to_x = np.asarray(angles, dtype='float')
            if (not self.angles_to_x_units) and np.any(self.angles_to_x):
                if np.max(self.angles_to_x) > 2 * np.pi:
                    print("\nSetting 'angles_to_x_units' to 'deg'")
                    self.angles_to_x_units = 'deg'
            self.intensities = np.asarray(intensities, dtype='float')
        else:
            print('Input file is not valid.')

    def save(self, fp, strict=False):
        """Save the data into a 'dfg' (diffractogram) file.

        Boolean 'strict' defines whether enforcing the file extension
        or printing a warning and continuing (default).

        """
        if type(fp) is str:
            fp = newname(fp)
            with open(fp, 'w') as f:
                self._write(f, strict)
                print('Data saved to file %s' % fp)
        elif hasattr(fp, 'closed'):
            # it's a file object
            if overwrite(fp.name):
                if fp.closed:
                    with open(fp.name, 'w') as f:
                        self._write(f, strict)
                        print('Data saved to file %s' % fp.name)
                else:
                    if 'w' in fp.mode:
                        self._write(fp, strict)
                        print('Data saved to file %s' % fp)
                    else:
                        raise IOError('%s is not open for overwriting' %
                                      fp.name)
    def _write(self, fp, strict):
        """Write data to OPEN file object 'fp'

        boolean 'strict' defines whether enforcing the file extension
        or printing a warning and continuing.

        """
        fpext = os.path.splitext(fp.name)[1]
        if fpext.lower() != Diffractogram.ext:
            msg = ('Warning. File extension should be %s instead of %s'
                   % (Diffractogram.ext, fpext))
            if strict:
                raise DeprecationWarning(msg)
            else:
                print(msg)
        headers = [
            '## structure: ' + self.structure,
            '## chem name: ' + self.chemname,
            '## space group: ' + self.sg_name + ' (' + str(self.sg_number) + ')',
            '## crystal system: ' + self.crystsys,
            ]
        ucst = '## unit cell: '
        if self.sep:
            sep = self.sep
        else:
            sep = ' '
        for i in self.unitcell:
            ucst += str(i) + sep
        ucst += (sep + '(' + self.uc_len_units + ')' +
                sep + '(' +  self.uc_ang_units + ')')
        headers.append(ucst)
        try:
            for line in headers:
                fp.write(line)
                fp.write(nl)
            fp.write(nl)
            col_head = ''
            cols = []
            if np.any(self.hkl):
                col_head += '##' + sep + 'h' + sep + '  k' + sep + '   l'
                cols.append(self.hkl)
            if np.any(self.d_spacings):
                col_head += sep + 'd_spacing' + sep + '(' + self.d_spacings_units +')'
                cols.append(self.d_spacings)
            if np.any(self.angles_to_x):
                col_head += sep + 'angle_to_x' + sep + '(' + self.angles_to_x_units +')'
                cols.append(self.angles_to_x)
            if np.any(self.intensities):
                col_head += sep + 'intensity'
                cols.append(self.intensities)
            fp.write(col_head)
            fp.write(nl)
            fp.write(nl)
            cols = np.asarray(cols).T
            for line in cols:
                fp.write(sep)
                line = str(line)       # numpy truncates with a good precision
                line = line.replace('(', '').replace(')', '')
                line = line.replace('[', '').replace(']', '')
                line = line.replace(',', '').replace('MillerPlane', '')
                line = line.split()     # make it a list
                for i in line[:3]:
                    # rigth justify hkl
                    i = i.rjust(3)
                    fp.write(i)
                    fp.write(sep)
                for i in line[3:]:
                    # right justify the rest
                    i = i.rjust(8)
                    fp.write(i)
                    fp.write(sep)
                fp.write(nl)
        finally:
            fp.close()

    def __repr__(self):
        out = self.__class__.__name__ + '{'
        if self.chemname:
            out += self.chemname + ' '
        if self.structure:
            out += '(' + self.structure + '), '
        if self.sg_name:
            out += self.sg_name
            if self.sg_number:
                out += ' (' + str(self.sg_number) + '), '
            else:
                out += ', '
        if self.crystsys:
            out += self.crystsys + ', '
        if np.any(self.unitcell):
            ucrepr = self.unitcell.__repr__().strip('array')
            out += ucrepr
        out += '}'
        if out == self.__class__.__name__ + '{}':
            out = self.__class__.__name__ + '{"' + self.filename + '"}'
        return out

    def __str__(self):
        out = 'Instance of <%s>\n\n' % self.__class__.__name__
        if self.filename:
            out += 'file name: %s\n' % self.filename
        if self.chemname:
            out += 'chemical name: %s\n' % self.chemname
        if self.structure:
            out += 'structure: %s\n' % self.structure
        if self.sg_name:
            out += 'space group: ' + self.sg_name
            if self.sg_number:
                out += ' (' + str(self.sg_number) + ')\n'
            else:
                out += '\n'
        if self.crystsys:
            out += 'crystal system: %s\n' % self.crystsys
        if np.any(self.unitcell):
            out+= 'unit cell: %s' % self.unitcell
        if self.uc_len_units and self.uc_ang_units:
            out += ' (%s) (%s)\n' % (self.uc_len_units, self.uc_ang_units)
        else:
            out += '\n'
        return out
